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We present density functional theory (DFT) calculations for 6if-SiC{0001} surfaces with differ- 
ent surface stackings and terminations. We compare the relative stability of different (0001) and 
(0001) surfaces in terms of their surface free energies. Removing surface and subsurface Si atoms, 
we simulate the formation of graphene and graphene-like overlayers by Si evaporation. We find 
that overlayers with a different nature of bonding are preferred at the two non-equivalent surface 
orientations. At (0001), a chemically bonded, highly strained and buckled film is predicted. At 
(0001), a van der Waals (vdW) bonded overlayer is preferred. We quantify the vdW binding and 
show that it can have a doping effect on electron behavior in the overlayer. 

PACS numbers: 68.65.Pq,73.22.Pr,73.20.At 



I. INTRODUCTION 

Graphene [1] is a highly interesting material for future 
nanoelectronic devices [2, 3] such as transistors [4], in- 
tegrated circuits [5] or detectors [6]. This is because of 
its unique electronic properties but also, and more im- 
portantly, because of the possibility to adjust and con- 
trol these properties. The nature and magnitude of con- 
duction and overall (opto-) electronic properties can be 
modified, for example, by applying electric fields [7], by 
adsorbants [8, 9], by utilizing finite-size effects [10, 11] or 
by an environment-induced material transformation into 
graphane [12, 13] or graphene oxide [14]. 

An important issue for the use of graphene and 
graphene-derivatives [15] for devices is the actual influ- 
ence of substrates (on which these materials are placed 
in the device) on the bandstructure and hence on elec- 
tron behavior. Although graphene and its derivatives do 
not tend to easily form chemical bonds to other mate- 
rials, van der Waals (vdW) interactions will always be 
present. For the fully hydrogenated graphene-derivative, 
graphane [12, 13], we [16] have recently reported first- 
principles vdW density functional (vdW-DF) calcula- 
tions [17], predicting that vdW interactions stabilize mul- 
tilayer formation. Moreover, we predict that the electron 
behavior in the graphane multilayer may deviate, at least 
locally, from the behavior in the monolayer [16]. Similar 
effects must naturally also be expected when graphene 
is simply exposed to a substrate. Indeed, combinations 
of vdW-DF and GW [18] calculations for graphene on 
various metal surfaces have already predicted that the 
vdW interaction can shift the graphene Fermi level [19] 
and that these shifts can be either positive or negative, 
depending on the the actual substrate material. 
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In this paper, we study graphene and graphene-like 
overlayers at 6H-SiC{0001} surfaces, focusing on vdW 
binding and the effect of vdW forces on electron behav- 
ior. The 6H-SiC{0001} surfaces are a natural choice for 
studying general substrate-graphene vdW interactions. 
SiC is regarded as promising substrate candidate for 
large-scale fabrication of pure graphene by Si evapora- 
tion from {0001} surfaces [20-22]. We first characterize 
the thermodynamic stability of SiC {0001} surfaces with 
different orientations, atomic stacking and surface ter- 
mination. We then simulate the Si evaporation be re- 
moving Si atoms from surface and subsurface layers, let- 
ting the systems find their new ground state. At (0001), 
we predict a chemically bonded, strongly buckled and 
stretched graphene-like overlayer. At (0001), we predict 
a flat vdW-bonded graphene overlayer. For the vdW- 
bonded overlayer we perform band-structure calculations 
and find a modified electron behavior indirectly induced 
through vdW forces. 

The paper is organized as follows. In Sec. II, we give 
a short background of SiC and its {0001} surfaces. We 
present details of our computational method in Sec. III. 
In Sec. IV we present our results for the stability of var- 
ious surfaces and surface/overlayer systems, vdW bind- 
ing between SiC and graphene and the resulting band 
structure. The results are discussed in Sec. V and we 
summarize and conclude our work in Sec. VI. 



II. MATERIALS BACKGROUND 

A. Bulk SiC 

Figure 1 details the atomic structure of 6i7-SiC. The 
crystal structure is hexagonal. Table I presents a compar- 
ison of the lattice parameters obtained from experiment 
[23] and from our first- principle calculations described 
further below. 



Along the c-axis, the bulk repeat unit is composed 
of six SiC bilayers. Each bilayer contains 50 % silicon 
and 50 % carbon and forms a buckled honeycomb lat- 
tice. We define the [0001] direction as indicated in the 
figure. With this definition, the Si atoms in each bilayer 
are located above (along [0001]) the center-of-mass plane, 
C atoms are located below. 

The stacking sequence (along [001]) is ABCA'C'B'. 
We use primes to distinguish the first A, 5, and C lay- 
ers in each repeat unit from the second ones. Because of 
the six layers in the repeat unit, naively one would ex- 
pect twelve possible surface terminations for each surface 
orientation. Due to symmetry, however, each two of the 
twelve terminations per surface are equivalent. A rota- 
tion by 180° around the c axis, maps A sites on A sites, 
B sites on C sites and C sites on B sites. Therefore, the 
rotated stacking is A'C B' ABC showing the equivalence 
of A and A^ sites, B and sites, and C and B^ sites. 



B. Surfaces 
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FIG. 1: (Colour online) Bulk structure of 6H-SiC. Si atoms 
are represented by yellow (light gray), large spheres; C atoms 
are represented by black (black), small spheres. The light 
gray dashed lines represent the unit cell. 



Figure 2 shows six out of the twelve different (ideal) 
6i^-SiC{0001} surface configurations. Along the [0001] 
direction, the 6H crystal structure lacks inversion sym- 
metry. Thus, the corresponding (0001) surface (the so- 
called nominally Si-terminated surface [24]) and (0001) 
surface (the so-called nominally C-terminated surface 
[24]) are different. For each surface orientation there 
are six different possible ideal terminations, correspond- 
ing to two chemically different terminations (Si or C) 
times three structurally different terminations. Ideal here 
means that we only consider full-coverage surfaces. In 
practice, a large surface may exhibit partial coverage to 
counteract a diverging surface dipole [25] and there may 
be surface reconstructions as to saturate surface dangling 
bonds. 

The top panels show the set of three most natural (fiat) 
Si-terminated (0001) surfaces. These are denoted by Sil, 
Si2, and Si3. In a Sil surface, the surface Si atoms are 
located on top of C atoms that are located in the first sub- 
surface SiC bilayer. In Si2 and Si3 surfaces, the surface 
Si atoms are correspondingly located on top of C atoms 
in the second and third subsurface SiC bilayer. The C- 
terminated counterparts (CI, C2, and C3) are obtained 
by adding an additional C layer on top of the surface Si 
layer. 

The set of bottom panels displays C-terminated (0001) 
surfaces. The three surfaces are denoted by CI, C2 and 





exp (Ref. 23) 


present 


a in A 


3.073 


3.091 


c in A 


15.118 


15.181 



TABLE I: Comparison of SiC bulk lattice parameters from 
experiment and our first-principle calculations. 



C3. In a CI surface, the surface C atoms are located on 
top of Si atoms that are located in the first subsurface 
SiC bilayer. In C2 and C3 surfaces, the surface C atoms 
are correspondingly located on top of Si atoms in the sec- 
ond and third subsurface SiC bilayer. The Si-terminated 
counterparts (Sil, Si2, and Si3) are obtained by adding 
an additional Si layer on top of the surface C layer. 

At the (0001) surface, C atoms bind only to a single 
Si atom below; they possess three dangling bonds. Si 
atoms, on the other hand, bind to three nearest-neighbor 
C atoms; they only possess one dangling bond. At (0001), 
the situation is reversed. Therefore, in the absence of 
complex reconstructions, the (0001) surface is intuitively 
expected to be Si terminated (and therefore denoted as 
the nominally Si-terminated surface [24]) and the (0001) 
surface is intuitively expected to be C terminated (and 
therefore denoted as the nominally C-terminated surface 
[24]). 



III. COMPUTATIONAL METHOD 
A. Surface and overlayer stability 

All surface and surface/overlayer calculations are per- 
formed with the planewave pseudopotential [26] DFT 
code dacapo [27] and the PBE [28] functional for ex- 
change and correlation. We use a planewave cutoff of 
400 eV and a (4 x 4 x 1) k-point sampling [29] Force 
relaxations are performed until the residual force is less 
than 0.03 eV/A. 

We use the supercell approach and represent each sur- 
face or surface/overlayer system by slab geometry. The 
supercells have a height of 40 A and the lateral dimen- 
sions are fixed to accommodate the (1x1) SiC {0001} 



3 



Sil 



Si2 



Si3 



T'T 1 r r 



r 




..jksBf '■'J^'^ ^^^^ ^ 



FIG. 2: (Color online) (1 x 1) 6H-SiC(0001) and (0001) sur- 
faces with different details stacking in the outermost layers. 
Color coding as in Fig. 1. The set of top panels shows the 
three nonequivalent (0001) surfaces with Si termination Sil, 
Si2 and Si3 (see text for a definition of this labeling). Corre- 
sponding C-terminated surfaces (CI, C2 and C3) are obtained 
by adding a C layer on top of terminating Si layer. The set 
of bottom panels shows the three nonequivalent (0001) sur- 
faces with C termination CI, C2 and C3. Corresponding Si- 
terminated surfaces (Sil, Si2 and Si3) are obtained by adding 
a Si layer on top of the terminating C layer. 



surface. Surfaces with different stackings and chemical 
terminations are represented by slabs of different thick- 
nesses, varying from 9 and 15 bilayers (plus one optional 
excess Si or C layer). Also, we saturate dangling bonds 
at the (0001) [(0001)] side of the slab that represents the 
(0001) [(0001)] surface by attaching H atoms. Because 
of the asymmetry of the slabs, we use a dipole correction 
[30]. With respect to thickness, the calculated surface- 
energy differences, see Table II, are converged by at least 
±2 meV/A^. We have tested this accuracy by compar- 
ing surface energies of equivalent surfaces that are repre- 
sented by slabs of different thicknesses. Specifically, any 
surface represented by a particular slab can be recovered 
by addition of three SiC bilayers. 



Thermodynamic comparison of chemically 
different terminations 



We determine the equilibrium surface preference by 
comparing surface free energies [31, 32]. Due to the lack 
of inversion symmetry along [0001], use of slab geometry 
only enables us to calculate the sum (denoted by a and 
defined below) of the two surface energies (denoted by 
a) corresponding to the (0001) surface and the (0001) 
surface. However, keeping the geometry and chemical 



composition fixed at one side of the slab, we ensure that 
the contribution to a from that side is always the same. 
We thus can compare the relative stability of different 
structures and compositions at the other side (not fixed) 
by considering the difference in a for various stackings 
and terminations. 

We define the sum of the two surface free energies as 
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^SiCSiC - (^C - ^Si)MC - riGT^GT . (1) 



Here, £^siab is the total energy of the 6i^-SiC surface slab, 
including a possible graphene-like overlayer, that con- 
tains nsi silicon atoms (per supercell), tiq carbon atoms 
that belong to the SiC and ncr carbon atoms that belong 
to the graphene overlayer. Furthermore, esic denotes the 
energy of one stoichiometric unit of bulk SiC, ecr is the 
energy of graphene,^ and /ic is carbon chemical poten- 
tial. We note that we have assumed equilibrium between 
the bulk and the surface, esic = Msi + Mc- 

The values of the carbon chemical potential in (1) are 
restricted to a finite range. If the carbon chemical po- 
tential is larger than the free energy per carbon atom in 
graphene (/ic > ^'graphene), the formation of (additional) 
graphene overlayers becomes more favorable. If, on the 
other hand, the silicon chemical potential is larger than 
the free energy per silicon atom in bulk Si (with diamond 
structure, jiisi > gsi)^ the formation of bulk Si is more 
favorable. Therefore, the allowed range of the carbon 
chemical potential is 



^graphene > > ^SiC " ^Si- 

We also introduce 
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(Jo, 



(2) 



(3) 



from which we infer the relative stability of surfaces with 
identical orientations but different stackings and chemi- 
cal compositions. 



C. van der Waals binding 

In Sec. IV B, we identify a chemically non-binding 
SiC (0001) /graphene system as thermodynamically sta- 
ble. We study this system further using the van der 
Waals density functional (vdW-DF) method [17]. In par- 
ticular we use the nonlocal correlation functionals E^^~^ 
and E^^~'^ of Refs. 33 and 34, respectively. 

We calculate the energy variation, including vdW 
forces, by a postprocessing method as follows. We first 
perform traditional DFT calculations (using the PBE 
exchange-correlation functional) for various separations 



^ This energy must contain the strain energy since the graphene 
overlayer is highly expanded in a SiC surface cell. 
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between the SiC and the graphene overlayer. In these 
calculations we choose the lateral dimensions of the unit 
cell to fit those of a (4 x 4) SiC(0001) surface. This allows 
us to study a (4 X 4) SiC(0001)/(5 x 5) graphene system 
in which the graphene overlayer is hardly strained at all. 
We use a (3 X 3 X 1) k-point sampling to ensure an accu- 
rate electronic density for further evaluations of nonlocal 
correlations according to 
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vdW-DF-?; 



(4) 



Here, E^^~'^[n] is the correlation energy from one of the 
non-local functionals of Refs. 33 {v = 1) and 34 [v = 2). 
Also, in Eq. (4), Eo_y[n] is given by 
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where E^"^^ is the VWN-LDA [35] correlation energy 

and El=^ = ^fvPBE ^v=2 _ ^rPW86 ^j^^ 

revPBE [36] and rPW86 (refitted form of PW86) [37]) 
exchange functionals. 

We determine the vdW binding separation and energy 
by finding the minimum in the layer-binding energy de- 
fined as 

^bind(G^) = ^vdW-DF(G^) " £^vdW-DF(G^ ^ Oo). (6) 

Here, d is the distance between the surface layer and 
the overlayer. For a detailed description of a robust im- 
plementation of the evaluation of Eq. (6), we refer to 
Refs. 38-41. 



D. Band-structure calculations 
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TABLE 11: Differences in surface energies for surfaces with 
identical orientations and terminations (which makes the 
numbers independent of the chemical potentials) but with dif- 
ferent detailed stacking of the outermost surface layers. The 
energetically most favorable surface for each orientation and 
termination is chosen as reference and therefore characterized 
by Act = 0. The effect of surface vibrations is not included. 
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1x1 (0001) 
73 X 73 (0001) 
3x3 (0001) 


-527 
-559 
-539 


-502 
-534 
-513 


-476 
-509 
-487 


1x1 (0001) 
\/3 X 73 (0001) 
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-128 
-75 
-118 


-102 
-49 
-93 


-77 
-23 
-67 



TABLE IIL Preference of surface chemical composition at 
6if-SiC{0001}. The table lists differences in surface free ener- 
gies between Si- and C-terminated surfaces for three values of 
the chemical potential of C. A negative value indicates higher 
stability of the Si-terminated surface. 



We also perform band- structure calculations for (5x5) 
graphene on (4 x 4) SiC(OOOl) to probe the effect of vdW 
bonding on electron behavior. The band-structure calcu- 
lations are performed as follows. We fix the SiC-graphene 
separation at the value predicted by the vdW-DF2 cal- 
culations and determine the density with a (4 x 4 x 1) 
k-point sampling to ensure a high accuracy in our large 
unit-cell calculations. This density is subsequently used 
to perform traditional GGA calculations of the energy 
spectra at various k points. We focus on the special 
Brillouin-zone points F = (0,0,0), K (2/3,1/3,0) and 
M = (1/2, 1/2,0) and the lines along KF, FM, and KM. 

Since we are interested in band-structure modifica- 
tions, we also perform the same type of calculations for a 
single graphene layer. In order to not encounter modifi- 
cations that may be solely due to a (very small) variation 
in the graphene lattice parameter, we use the same unit 
cell as in the case for the SiC/graphene system. 

IV. RESULTS 

A. Relative stability of clean 6i/-SiC{0001} surfaces 

a. Stacking preference for fixed chemical termination. 
Table II compares surface energies for surfaces with 



identical orientations and identical chemical composi- 
tions but with a different detailed stacking of the surface 
layers. For each orientation and chemical composition of 
the surface, we choose the configuration with lowest sur- 
face energy, or more precisely with lowest cr, as reference. 
This configuration is thus characterized by Act = 0. 

We find that the surface-energy differences are of the 
order of a few meV/A^. These values are close to or 
below the actual accuracy of our calculations. Thus, the 
surface-energy differences between surfaces with different 
stackings are too small to be of significant importance for 
determining the stable surface configuration. 

h. Chemical composition. In Table III we deter- 
mine the preferred chemical composition of the (0001) 
and (0001) surfaces. For each orientation we list dif- 
ferences in surface free energies between Si- and C- 
terminated surfaces at three different values of the C 
chemical potential. The three values of the C chemical 
potential correspond to its maximal value, its minimal 
value, and the value in between. 

We find that, within the entire allowed range of the C 
chemical potential, the Si-terminated surfaces possesses 
lowest free surface energy, independent of the orientation. 
For (0001) orientation, this coincides with the expecta- 
tions on the basis of the number of dangling bonds. For 
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FIG. 3: (Color online) Formation of graphitic overlayers at SiC{0001}. Color coding as in Fig. 1; in addition, C atoms forming 
an overlayer are in red (dark gray) and connected through bonds. The set of top panels shows side views on truncated (0001) 
and (0001) surface structures with either Si or C termination. The set of mid-panels shows side views the same structures 
(i.e. unrelaxed), but with (sub) surface Si layers removed. The set of bottom panels shows side and top views on the 
geometries obtained by relaxing the systems in the mid-panels. At each surface there exist chemically bonded systems (first 
and third column in the bottom set of panels) and systems that are not chemically bonded (second and fourth set of panels). 
Our thermodynamics analysis indicates that the chemically bonded system is the stable one at SiC(OOOl); At SiC(OOOl), the 
chemically non-binding configuration is stabilized. The stable, that is, thermodynamically favored configurations are underlined 
in green (gray). 
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(0001) surfaces, the predictions appear counterintuitive. 
There, the Si-terminated surface possesses more danghng 
bonds than the C-terminated surface and should there- 
fore be less favorable. 

Table III also reports that the perhaps counterintu- 
itive result for (0001) surfaces remains unaltered when 
considering larger surface unit-cells that allow for sur- 
face reconstructions (without considering more complex 
surface terminations than pure Si or C termination, how- 
ever), such as (\/3 X ^/3) and (3 x 3) surface unit-cells. 
Details of our calculations concerning surface reconstruc- 
tions are documented in the supplementary material. 

The appearant equilibrium preference of Si-terminated 
(0001) surfaces reflects the rich phase diagram of SiC. 
Our comparison, limited to full Si and C coverages, is 
likely not exhaustive enough to capture this richness in a 
more quantitative manner. On the other hand, resolving 
the difficult structure is not relevant for our search for 
carbon overlayers at SiC {0001} surfaces that may form 
by evaporation of Si atoms. 



B. Graphitic overlayer formation by Si evaporation 

We now turn our focus towards structure, stability and 
bonding of graphitic overlayers at SiC{0001}. Figure 3 
illustrates the formation of different graphitic overlayers 
at SiC {0001} surfaces by evaporation of Si atoms [20-22]. 
We use (1x1) surface unit-cells to study the effect of re- 
moving full surface and subsurface layers. As illustrated 
in the preceding subsection and the supplementary ma- 
terial, the initial surface morphologies can, in principle, 
be much more complicated. However, since the surface 
atoms evaporate in a heating process, we expect that the 
detailed surface structure before heating is of minor im- 
portance. 

a. Structure. The set of top panels of Fig. 3 shows 
four (1x1) surfaces with different orientations and termi- 
nations. In the set of mid-panels, we have removed sur- 
face and subsurface Si atoms. The set of bottom panels 
presents the geometries that are obtained by relaxing the 
structures in the set of mid-panels above. The figure also 
illustrates that the role of the surface composition (and 
structure) effectively reduces to determining the amount 
of Si that needs to be evaporated to generate a specific 
final SiC/graphene system. 

At both surface orientations we obtain two types of 
overlayers. The first type of overlayer is chemically 
bonded, see ffist and third column in the set of bottom 
panels in Fig. 3. As a result the overlayer adopts the SiC 
lattice parameter resulting in a considerable stretching 
of the C bonds. At (0001), the overlayer preserves the 
hexagonal graphene-like shape but also exhibits a buck- 
ling. At (0001), the hexagonal shape is not preserved. 
The C atoms arrange themselves in chains located in 
bridge positions at the surface. 

The second type of overlayer is not chemically bonded, 
see second and fourth column in the set of bottom panels 
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FIG. 4: Energy variation of graphene on SiC (0001) including 
vdW forces. 



of Fig. 3. The non-binding character is reflected by the 
large separations between the overlayers and the outer- 
most surface layers. In both cases, the absence of chem- 
ical bonding leaves the ideal hexagonal graphene shape 
unchanged. We expect that these carbon overlayers con- 
tract and adopt the unstrained graphene lattice. 

b. Stability. We compare the stability of the differ- 
ent surface/overlayer systems with identical orientations 
by means of their surface energy, see Eq. (1). For the 
chemically bonded systems, we consider the overlayer to 
be in equilibrium with the SiC surface, that is, the num- 
ber of graphene units is set to zero. For the other two 
systems, we assume that the chemical potential of the 
overlayer is not related to that of SiC, but to that of 
(strained) graphene (so correcting for the strain energy 
due to the lattice misfit). 

At the (0001) surface, we find that the chemically 
bonded system is more preferred. The surface-energy dif- 
ferences are 86 meV/A^ at the minimal value of carbon 
chemical potential and 137 meV/A^ at the maximal value 
of carbon chemical potential. At the (0001) surface, the 
chemically non-bonding system is preferred. The corre- 
sponding surface-energy differences are 445 meV/A^ (at 
/i^^^) and 393 meV/A^ (at /i^^^). 



C. vdW bonded (5 x 5) graphene at (4 x 4) SiC(OOOl) 

We quantify the vdW binding of graphene at SiC (000 1) 
by studying the vdW energy variation of the system 
shown in the rightmost bottom panel of Fig. 3. In the 
depicted (1x1) system, the graphene overlayer is highly 
strained. We therefore use a (4 x 4) SiC(OOOl) surface 
unit- cell in our calculations combined with a (5x5) 
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FIG. 5: Zone folding of the graphene band diagram along 
FM. The left panel shows the unfolded calculated (1x1) 
band structure of graphene. The diagram in the mid-panel is 
constructed from that in the left by folding five times. The 
right panel shows the calculated band diagram using a (1 x 5) 
graphene unit cell. 



graphene overlayer [42] to ensure that the overlayer is 
almost unstrained (the C-C bond length is stretched by 
- 0.5% only). 

Figure 4 shows the energy variation as a function of 
the separation between the SiC surface-layer and the 
graphene overlayer. Both, results using vdW-DFl and 
vdW-DF2 are shown. The insert shows that traditional 
GGA calculations (PBE) do not predict any meaningful 
binding (notice the different scale on the y-ax.is in the 
insert). 

The vdW-DF energy variations agree qualitatively. 
vdW-DF2 predicts a slightly smaller binding separation 
and a slightly larger binding energy. The numerical val- 
ues are (ibind = 3.6 A and £^bind = —54.7 meV per carbon 
atom (in graphene) for vdW-DFl and (ibind = 3.4 A and 
^bind = —55.6 meV per carbon atom for vdW-DF2. 

We note that accounting for vdW binding for the chem- 
ically non-binding graphene overlayer at SiC (0001) (sec- 
ond panel from the left in bottommost set of panels in 
Fig. 3) is expected to lower surface energy by a similar 
amount. However, our calculated values correspond to 
^vdw^ ^ 13 meV/A. Therefore, even with an account of 
vdW interactions, the chemisorbed carbon overlayer will 
still be more preferable than the chemically nonbinding 
(vdW-bonded) overlayer. 



the band structure of graphene (without substrate) along 
a straight line from F to K as calculated within a (1 x 
1) unit cell. In the mid-panel we show a band diagram 
that is constructed from the left panel by zone folding it 
five times. The right panel shows the band diagram as 
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D. Band structure of vdW-bonded graphene at 
nominally C-terminated SiC (0001) 

a. Consistency check for zone folding. In Fig. 5 we 
check that our large unit-cell calculations capture and 
reliably reproduce the details of the electron behavior, 
that is, the band-structure physics. The left panel shows 



FIG. 6: (Color online) Electronic structure of vdW bonded 
(5 X 5) graphene at (4 x 4)SiC(000i). The top panel shows 
the calculated band diagram along KF, FM and MK. The 
UVB and LCB corresponding to graphene are highlighted. 
The set of bottom panels shows (absolute values of) various 
wave functions, (a) Graphene UVB WF in K, (b) graphene 
LCB WF in K, (c) WF in ki (see top panel) and (d) WF in 
(see top panel). 
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calculated within a (1 x 5) unit cell. 

The constructed and the calculated zone-folded band 
diagram agree reasonably well. The slight differences 
may be due to differences in the underlying electronic 
densities that are used in the respective calculations and 
which is transferred from the (1x1) band diagram to the 
zone-folded diagram. 

b. Overlayer band- structure. The top panel of 
Fig. 6 shows the calculated band diagram of (5 x 5) 
graphene at (4 x 4) SiC(OOOl). We restrict the plot to 
relevant energy window around the Fermi level. The com- 
plex band structure is due to the (5 x 5) zone folding of 
graphene bands and the (4 x 4) zone folding of SiC bands. 

Among the many bands, two bands are highlighted. 
These bands correspond to the upper valence band 
(UVB) and to the lower conduction band (LCB) in an 
isolated graphene sheet. We have explicitly checked that 
the wave functions (WF) corresponding to the UVB and 
LCB are localized on the graphene overlayer. 

The set of bottom panels of Fig. 6 shows various WFs. 
These WFs are representative for the different types of 
WF localization in the system. Panels (a) and (b) show 
that the UVB and LCB WFs in K, for example, are en- 
tirely located on the overlayer. WFs fully localized on 
graphene are typical for the LCB. 

We note that a kink arises in the band structure vari- 
ation as the graphene band crosses the Fermi level at ki 
(see top panel of Fig. 6 for a definition of ki). There 
we find that the graphene UVB WF can also be shared 
between the SiC and the graphene. 

Finally, in panel (d), we illustrate that the graphene 
overlayer also slightly affects the nature of the SiC states 
at the Fermi level in k2 (see top panel of Fig. 6 for a 
definition of k2)- Although most WFs corresponding to 
a band between the UVB and LCB are fully localized 
within the SiC substrate, at kd and other k points, some 
WFs do have a small weight also on the graphene over- 
layer. 



V. DISCUSSION 
A. Surface stability 

Our calculations of (1 x 1) SiC{0001}, including some 
(v^ X \/3) and (3 x 3) reconstructions in the supplemen- 
tary material, predict a preference for Si termination. 
This prediction is reasonable for the (0001) surface but 
surprising for the (0001) surface for which we would ex- 
pect a C termination. 

Resolving this discrepancy requires a more careful 
study that take into account the actual growth condi- 
tions of SiC [43] and/or improves on the description of 
surface reconstructions. The latter task would have to 
consider a richer set of surface terminations also includ- 
ing excess and deficiency Si or C. However, the size of 
the {^/3 X \/3) and (3 x 3) surface unit-cells makes a full 



reconstruction-search a large project of its own, requir- 
ing systematic structure-search strategies such as consid- 
ering a larger pool of candidate geometries with struc- 
tural motifs [44] of the reconstructed surfaces obtained 
here, evolutionary- type of iteration [45] or other global 
structure-search methods [46-48]. 

Such a search for surface reconstructions is clearly be- 
yond the scope of the present work, in particular, since 
the exact morphology of the stable SiC {0001} surfaces 
is only of minor importance for the main objective of 
this paper: the study of modification of electron behav- 
ior in graphene over layers due to vdW interactions. At 
the same time, we emphasize that our calculated surface 
energies are upper limits of the true surface energies. 



B. Nature of binding in SiC/graphene systems 

We have identified preferred SiC/graphene systems 
as they may result by Si evaporation from various 
SiC{0001} surfaces. Our results indicate that the nature 
of binding at the nominally Si-terminated (0001) surfaces 
is different from the nature of binding at the nominally 
C-terminated (0001) surface. 

At (0001), we have identified a chemisorbed, strongly 
buckled graphene overlayer. At (0001), the overlayer is 
stabilized by vdW forces. The different nature of binding 
at the two surfaces may have consequences for the qual- 
ity of graphene that is grown by evaporation of Si from 
different SiC {0001} faces (at the nominally C-terminated 
or at the nominally Si-terminated face). 



C. Electron behavior in vdW-bonded graphene 

Our band-structure calculations for SiC/graphene, see 
Fig. 6, show that vdW binding can cause a doping-like 
effect. In free-standing graphene, the density of states 
(DOS) vanishes at the Fermi level. ^ On SiC, the vdW 
binding renders graphene a p-doped metal. The Fermi 
level is shifted to lower energies into the original valence 
band where the DOS is finite. The prediction of a Fermi- 
level shift is similar to the results reported for graphene 
on various metal surfaces in Ref. 19. 

In the present case, vdW binding leads to a further 
modification in the band-structure. At the new Fermi 
level between K and F, we also observe a kink in the band 
structure. To the left of ki (see Fig. 6) the dispersion 
in the LCB can be fitted to a parabolic form E{k) = 
E'o + a- (A: — /ci)^, from which we infer an effective electron 



We have checked that this remains true also if we shghtly stretch 
the C-C bond length such that a (5 X 5) unit cell is commensurate 
with a (4 X 4) SiC{0001} cell. 
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mass 



meff = h 



(7) 



meff ~ 1.02 X 10 kg or meff 
(actual) electron mass. 



1.1 me where is the 



VI. SUMMARY AND CONCLUSIONS 

We present DFT calculations for SiC{0001} surfaces 
and graphitic overlayers at SiCjOOOl} surfaces. In partic- 
ular, we focus on surface and surface/overlayer stability, 
and on binding and band-structure modifications (due to 
the binding) of the overlayers. 

For surfaces, we study the relative stability as func- 
tion of the detailed stacking, as function of the chemical 
composition and to some extent (see supplementary ma- 
terial), as function of the type of reconstruction. We find 
that the surface-energy differences due to different de- 
tailed stacking of the outermost surfaces are below the 
accuracy of our calculations and not of any significance. 
For ideally truncated surfaces (no excess or deficiency of 
Si or C) we find that Si-terminated surfaces are generally 
more favorable than C- terminated surfaces. 



For SiC/overlayer systems we find two different types 
of overlayers. At SiC(OOOl), we predict an overlayer 
that is chemically bonded to the substrate. Because of 
the chemical bonding this overlayer is expected to sig- 
nificantly differ in its electronic nature from single-layer 
graphene. At SiC(OOOl), we predict a vdW-bonded over- 
layer. 

In line with Ref. 19, our band- structure calculations 
for the vdW-bonded graphene show that vdW interac- 
tions with a substrate can have a doping effect (here: p 
doping). As a novel feature, we also identify a kink in 
the electron dispersion at the Fermi level and calculate 
an effective mass of meff ^ l.lme at the minimum of the 
conduction band at this kink. 
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Relative stability of 6i7-SiC{0001} surface 
terminations and formation of graphene overlayers 
by Si evaporation: Supplementary material 

In the main article entitled 'Relative stability of 
6i^-SiC{0001} surface terminations and formation of 
graphene overlayers by Si evaporation', we investigate, 
among others, the relative stability of 6i^-SiC{0001} sur- 
faces as a function of their surface termination. Focusing 
on (1 X 1) surface unit-cells, we find that, independent 
of the orientation. Si termination is preferred over C ter- 
mination. This result seems to be partially in confiict 
with intuition based on counting the number of dangling 
bonds; at SiC(OOOl), a C-terminated surface would be 
expected. 

Here we present calculations where we consider the 
larger {^/3x^/3) and (3x3) surface unit-cells. These cells, 
in principle, allow for various surface reconstructions, and 
the absence of such reconstructions in (1x1) surface unit- 
cells could be the origin of the counterintuitive prediction 
of a preferred Si termination at SiC(OOOl). 

For all calculations we use a common 400 eV planewave 
cutoff. For {^/3 x ^/3) surface cells, a (3 x 3 x 1) k-point 
sampling is used; for (3 x 3) surface cells, a (2 x 2 x 1) k- 
point sampling is used. The surfaces are modeled as slabs 
consisting of at least six SiC bilayers. All slabs are ideally 
truncated, that is, they possess either a full-coverage Si- 
terminated surface or a full-coverage C-terminated sur- 
face. The surface models are relaxed until the forces on 
atoms no longer exceeds 0.01 eV/A. 

Figure 7 details different (3 x 3) surface reconstruc- 
tions obtained from our calculations. In fact, the direct 
relaxation of truncated surface slabs only produces fiat 
surfaces without true reconstructions, the Si-terminated 
(3 X 3) (0001) surface being an exception. All other re- 
constructions are obtained by pulling one atom out of 
the fiat surface-layer by 0.5 A and then restarting the 
relaxations. 

The most pronounced reconstruction, giving rise to tri- 
angular features, is identified at the C-terminated (3 x 3) 
(0001) surface, see set of mid-panels in Fig. 7. At Si- 
terminated (3 X 3) (0001), see set of top panels in Fig. 7, 
and Si-terminated (3 x 3) (0001), see set of bottom pan- 
els in Fig. 7, only smaller departures from the fiat surface 
are found. The C-terminated (3 x 3)SiC(000l), remains 
unreconstructed even after triggering reconstructions by 
pulling atoms out of the surface. 

For (a/3 X a/3) surface unit-cells, all but the Si- 
terminated (0001) configuration remain fiat. The recon- 
struction of the Si-terminated (a/3 x a/3) (0001) surface is 
similar to the reconstruction of the Si-terminated (3 x 3) 
(0001) surface and therefore not shown. 

The reconstructed configurations generally possess 
lower energies than their unreconstructed, fiat counter- 
parts. These energies are used in our comparison of sur- 
face energies in Table III in the paper. 

We are aware of the likely fact that the here-presented 
reconstructions only represent a small sample of the rich 
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FIG. 7: Top and side views of (3 x 3) surface reconstructions 
identified in this study. Different colors are used to distinguish 
surface atoms that undergo strong rearrangements from those 
in subsurface layers. 
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SiC surface phase-diagram. The {V3 X V3) and (3 x 3) 
surface unit-cehs ahow for a large freedom in the varia- 
tion of the coverage of surface Si or C. Closely related 
is the problem of identifying the lowest-energy structure 
for each coverage. A full reconstruction search would re- 
quire more advanced strategies and methods than used 
here and is outside the scope of the present work. 

Finally, we notice that identification of the full spec- 



trum of surface reconstructions is in fact only of sec- 
ondary relevance for the actual purpose of the main pa- 
per: the study of band-structure modifications due to van 
der Waals binding of graphitic overlayers at SiC{0001}. 
These overlayers arise from evaporation of Si atoms from 
SiC{0001}. As discussed in the main paper, the detailed 
structure and composition of the surface prior to evapo- 
ration may be only of minor relevance. 



